﻿



library(PBSmodelling)
library(BRugs)
library(R2OpenBUGS)
library(coda)


coda.data<-read.openbugs(stem="")
results<-as.data.frame(rbind(coda.data[[1]], coda.data[[2]], coda.data[[3]]))



# step 1: draw 1000 sets of paramaters
# step 2a: retrieve vector of 1000 conditional values for each X
# step 2b: take 25th highest and lowest lines at each of 1000 points along x axis
# step 3: Plot: draw lines through these points, shade in area, draw in expected value curve



x<-(seq(-1,1,length.out=1000))
y<-(seq(-1,1,length.out=1000))

## NOTE: the data used in the original model are coded so liberal is "high" and conservative "low," so that puts liberals on the right.
##       This main analysis file has lines that indicate a "flip" in the comment -- I am just rotating the figure so that liberal is on the left.
##       This translates to a full 180 degree rotation in both the x axis (about itself) and then the y axis. 

# liberals
upper.l<-rep(NA, 1000)
lower.l<-rep(NA, 1000)
mu.bar<-(-1.002)
for (j in 1:1000) {
mu.bar<-(mu.bar+0.002)
r<-rep(NA, 1000)
for (i in 1:1000) {
random<-as.integer(runif(1,1,length(results$a.1)))  # select a draw from posterior
curve<-function(x) {
results$rho.1[random]*x+results$rho.4[random]*x^2+results$rho.7[random]+results$rho.8[random]+results$rho.5[random]*x^2}
r[i]<-curve(mu.bar)
}
upper.l[j]<-sort(r)[975]
lower.l[j]<-sort(r)[25]
}
upper.l.curve<-lowess(seq(-1,1, length.out=1000), upper.l, f=.1)
lower.l.curve<-lowess(seq(-1,1, length.out=1000), lower.l, f=.1)

## flip the y-axis, first reversing the values
upper.l.curve$y<-rev(upper.l.curve$y)
lower.l.curve$y<-rev(lower.l.curve$y)
## then correcting the sign
upper.l.curve$y<-(-1)*upper.l.curve$y
lower.l.curve$y<-(-1)*lower.l.curve$y


# moderates
upper.m<-rep(NA, 1000)
lower.m<-rep(NA, 1000)
mu.bar<-(-1.002)
for (j in 1:1000) {
mu.bar<-(mu.bar+0.002)
r<-rep(NA, 1000)
for (i in 1:1000) {
random<-as.integer(runif(1,1,length(results$a.1)))  # select a draw from posterior
curve<-function(x) {
results$rho.1[random]*x+results$rho.7[random]+results$rho.5[random]*x^2}
r[i]<-curve(mu.bar)
}
upper.m[j]<-sort(r)[975]
lower.m[j]<-sort(r)[25]
}
upper.m.curve<-lowess(seq(-1,1, length.out=1000), upper.m, f=.1)
lower.m.curve<-lowess(seq(-1,1, length.out=1000), lower.m, f=.1)

## flip the y-axis, first reversing the order of the values
upper.m.curve$y<-rev(upper.m.curve$y)
lower.m.curve$y<-rev(lower.m.curve$y)
## then correcting the sign
upper.m.curve$y<-(-1)*upper.m.curve$y
lower.m.curve$y<-(-1)*lower.m.curve$y


# conservatives
upper.c<-rep(NA, 1000)
lower.c<-rep(NA, 1000)
mu.bar<-(-1.002)
for (j in 1:1000) {
mu.bar<-(mu.bar+0.002)
r<-rep(NA, 1000)
for (i in 1:1000) {
random<-as.integer(runif(1,1,length(results$a.1)))  # select a draw from posterior
curve<-function(x) {
results$rho.1[random]*x+results$rho.6[random]*x^2+results$rho.7[random]+results$rho.9[random]+results$rho.3[random]+results$rho.5[random]*x^2}
r[i]<-curve(mu.bar)
}
upper.c[j]<-sort(r)[975]
lower.c[j]<-sort(r)[25]
}
upper.c.curve<-lowess(seq(-1,1, length.out=1000), upper.c, f=.1)
lower.c.curve<-lowess(seq(-1,1, length.out=1000), lower.c, f=.1)

## flip the y-axis, first reversing the order of the values
upper.c.curve$y<-rev(upper.c.curve$y)
lower.c.curve$y<-rev(lower.c.curve$y)
## then correcting the sign
upper.c.curve$y<-(-1)*upper.c.curve$y
lower.c.curve$y<-(-1)*lower.c.curve$y


# FIGURE 2
op<-par(mfrow=c(1,3), mar=c(4,4,2,1)+0.1, oma=c(0,0,4,0))
plot(x,y, type = "n", main=paste("Liberals"), xlab=("(- Lib.)       Table Mean Ideology       (+ Cons.)"), ylab=("(- Liberal)          ∆θ         (+ Conservative)"))# setting up coord. system
polygon(c(x,rev(x)), c(lower.l.curve$y, rev(upper.l.curve$y)), col="lightcyan", border=NA)
lines(upper.l.curve, type="l", col="lightsteelblue", lty=2)
lines(lower.l.curve, type="l", col="lightsteelblue", lty=2)
curve<-function(x) {
mean(results$rho.1)*x+mean(results$rho.4)*x^2+mean(results$rho.7)+mean(results$rho.8)+mean(results$rho.5)*x^2}
curve1<-function(x) {
mean(results$rho.1)*x+(-2)*mean(results$rho.4)*x^2+mean(results$rho.7)+mean(results$rho.8)+mean(results$rho.5)*x^2}
mu.x<-curve(seq(-1,1,length.out=1000))
## flip correctly
mu.x<-(-1)*rev(mu.x)
mu.x1<-curve1(seq(0,1,length.out=1000))
## flip correctly
mu.x1<-(-1)*rev(mu.x1)
lines((seq(-1,1,length.out=1000)), mu.x, type="l", lwd=3, col="steelblue")
lines((seq(-1,0,length.out=1000)), mu.x1, type="l", lwd=3, col="steelblue", lty=4)
text(0.2,-0.8, "Hypothetical")
arrows(-0.1,-0.8,-0.65,-0.9, length=0.07, angle=20)
text(-0.4,0.3,"Observed")
arrows(-0.4,0.25,-0.78,-0.41, length=0.07, angle=20)

plot(x,y, type = "n", main=paste("Moderates"), xlab=("(- Lib.)       Table Mean Ideology       (+ Cons.)"), ylab=("(- Liberal)          ∆θ         (+ Conservative)"))# setting up coord. system
polygon(c(x,rev(x)), c(lower.m.curve$y, rev(upper.m.curve$y)), col="grey90", border=NA)
lines(upper.m.curve, type="l", col="grey70", lty=2)
lines(lower.m.curve, type="l", col="grey70", lty=2)
curve<-function(x) {
mean(results$rho.1)*x+mean(results$rho.7)+mean(results$rho.5)*x^2}
mu.x<-curve(seq(-1,1,length.out=1000))
## flip correctly
mu.x<-(-1)*rev(mu.x)
lines((seq(-1,1,length.out=1000)), mu.x, type="l", lwd=3, col="grey30")

plot(x,y, type = "n", main=paste("Conservatives"), xlab=("(- Lib.)       Table Mean Ideology       (+ Cons.)"), ylab=("(- Liberal)          ∆θ         (+ Conservative)"))# setting up coord. system
polygon(c(x,rev(x)), c(lower.c.curve$y, rev(upper.c.curve$y)), col="mistyrose", border=NA)
lines(upper.c.curve, type="l", col="mistyrose3", lty=2)
lines(lower.c.curve, type="l", col="mistyrose3", lty=2)
curve<-function(x) {
mean(results$rho.1)*x+mean(results$rho.6)*x^2+mean(results$rho.7)+mean(results$rho.9)+mean(results$rho.3)+mean(results$rho.5)*x^2}
curve1<-function(x) {
mean(results$rho.1)*x+(-2)*mean(results$rho.6)*x^2+mean(results$rho.7)+mean(results$rho.9)+mean(results$rho.3)+mean(results$rho.5)*x^2}
mu.x<-curve(seq(-1,1,length.out=1000))
## flip correctly
mu.x<-(-1)*rev(mu.x)
mu.x1<-curve1(seq(-1,0,length.out=1000))
## flip correctly
mu.x1<-(-1)*rev(mu.x1)
lines((seq(-1,1,length.out=1000)), mu.x, type="l", lwd=3, col="maroon")
## change the domain correctly for the hypothetical line, for flipped it would be 0,1
lines((seq(0,1,length.out=1000)), mu.x1, type="l", lwd=3, col="maroon", lty=4)
text(-0.2,0.8, "Hypothetical")
arrows(0.09,0.8,0.8,0.9, length=0.07, angle=20)
text(0.4,-0.45,"Observed")
arrows(0.4,-0.4,0.7,0.35, length=0.07, angle=20)
title("Ideological Persuasion: No Evidence of Polarization", line=2, outer=TRUE, cex.main=2)
par(op)



# FIGURE 3



# liberals
upper.l<-rep(NA, 1000)
lower.l<-rep(NA, 1000)
mu.bar<-(0.5)
for (j in 1:1000) {
mu.bar<-(mu.bar+0.002)
r<-rep(NA, 1000)
for (i in 1:1000) {
random<-as.integer(runif(1,1,length(results$a.1)))  # select a draw from posterior
curve<-function(x) {
-1*(results$rho.7[random]*x+results$rho.8[random]*x +results$rho.2[random])}
r[i]<-curve(mu.bar)
}
upper.l[j]<-sort(r)[975]
lower.l[j]<-sort(r)[25]
}
upper.l.curve<-lowess(seq(0.5,2.5, length.out=1000), upper.l, f=.1)
lower.l.curve<-lowess(seq(0.5,2.5, length.out=1000), lower.l, f=.1)





# moderates
upper.m<-rep(NA, 1000)
lower.m<-rep(NA, 1000)
mu.bar<-(0.5)
for (j in 1:1000) {
mu.bar<-(mu.bar+0.002)
r<-rep(NA, 1000)
for (i in 1:1000) {
random<-as.integer(runif(1,1,length(results$a.1)))  # select a draw from posterior
curve<-function(x) {
-1*(results$rho.7[random]*x)}
r[i]<-curve(mu.bar)
}
upper.m[j]<-sort(r)[975]
lower.m[j]<-sort(r)[25]
}
upper.m.curve<-lowess(seq(0.5,2.5, length.out=1000), upper.m, f=.1)
lower.m.curve<-lowess(seq(0.5,2.5, length.out=1000), lower.m, f=.1)



# conservatives
upper.c<-rep(NA, 1000)
lower.c<-rep(NA, 1000)
mu.bar<-(0.5)
for (j in 1:1000) {
mu.bar<-(mu.bar+0.002)
r<-rep(NA, 1000)
for (i in 1:1000) {
random<-as.integer(runif(1,1,length(results$a.1)))  # select a draw from posterior
curve<-function(x) {
-1*(results$rho.7[random]*x+results$rho.9[random]*x +results$rho.3[random])}
r[i]<-curve(mu.bar)
}
upper.c[j]<-sort(r)[975]
lower.c[j]<-sort(r)[25]
}
upper.c.curve<-lowess(seq(0.5,2.5, length.out=1000), upper.c, f=.1)
lower.c.curve<-lowess(seq(0.5,2.5, length.out=1000), lower.c, f=.1)



# SD effect: 
x<-(seq(0.5,2.5,length.out=1000))
y<-(seq(-1.5,1.5,length.out=1000))

# Note that multiplying by (-1) below does the correct flip if the policy indicators are NOT reverse coded
plot(x,y, type = "n", main=paste("Effect of Disagreement on Preference Direction"), xlab=("Table SD Ideology"), ylab=("(- Lib.)          ∆θ         (+ Cons.)"))# setting up coord. system
polygon(c(x,rev(x)), c(lower.l.curve$y, rev(lower.c.curve$y)), col="lightcyan", border=NA)
lines(upper.l.curve, type="l", col="lightsteelblue", lty=2)
lines(lower.l.curve, type="l", col="lightsteelblue", lty=2)
polygon(c(x,rev(x)), c(upper.l.curve$y, rev(upper.c.curve$y)), col="mistyrose", border=NA)
lines(upper.c.curve, type="l", col="mistyrose3", lty=2)
lines(lower.c.curve, type="l", col="mistyrose3", lty=2)
polygon(c(x,rev(x)), c(lower.c.curve$y, rev(upper.l.curve$y)), density=25, col="mediumpurple1", border=NA)
curve<-function(x) {-1*(
mean(results$rho.7)*x+mean(results$rho.8)*x+ mean(results$rho.2))}
mu.x<-curve(seq(0.5,2.5,length.out=1000)) 
lines((seq(0.5,2.5,length.out=1000)), mu.x, type="l", lwd=3, col="steelblue")
curve<-function(x) {-1*(
mean(results$rho.7)*x+mean(results$rho.9)*x+mean(results$rho.3))}
mu.x<-curve(seq(0.5,2.5,length.out=1000)) 
lines((seq(0.5,2.5,length.out=1000)), mu.x, type="l", lwd=3, col="maroon", lty=5)
curve<-function(x) {-1*(
mean(results$rho.7)*x+0)}
mu.x<-curve(seq(0.5,2.5,length.out=1000)) 
lines((seq(0.5,2.5,length.out=1000)), mu.x, type="l", lwd=3, col="grey30", lty=4)
text(x=1.3,y=1.3, "Liberal/Conservative Confidence Interval Overlap")
arrows(1.4,1.23,2.3,-0.15, length=0.07, angle=20)
legend(x=0.5, y=-1.0, bty="n", legend=paste(c("Conservative", "Moderate", "Liberal")), fill=c("maroon", "grey30", "steelblue"))


# TEST: liberal and conservative slopes in expected direction but not statistically different from each other....
plot(density(results$rho.7+results$rho.8, kernel="gaussian"), main=paste("Liberals"), xlab=(""), ylab=("Density"), xlim=c(-0.5,0.5))
lines(density(results$rho.7+results$rho.9, kernel="gaussian"), main=paste("Liberals"), xlab=(""), ylab=("Density"), xlim=c(-0.5,0.5), col="maroon")




# 




# FIGURE 4
### DELIBERATIVE PERSUASION GRAPH
op<-par(mfrow=c(6,3), mar=c(4,4,2,1)+0.1, oma=c(0,0,5,0))
plot(density(results$s.1+results$s.2, kernel="gaussian"), main=paste("Liberals"), xlab=(""), ylab=("Cut Programs"), xlim=c(-0.25,1))
abline(v=0, col="grey", lty=3)
abline(v=mean(results$s.1+results$s.2), col="black", lty=1)
plot(density(results$s.2, kernel="gaussian"), main=paste("Moderates"), xlab=(""), ylab=(""), xlim=c(-0.25,1))
abline(v=0, col="grey", lty=3)
abline(v=mean(results$s.2), col="black", lty=1)
plot(density(results$s.3+results$s.2, kernel="gaussian"), main=paste("Conservatives"), xlab=(""), ylab=(""), xlim=c(-0.25,1))
abline(v=0, col="grey", lty=3)
abline(v=mean(results$s.3+results$s.2), col="black", lty=1)

plot(density(results$d.1+results$d.2, kernel="gaussian"), main=paste(""), xlab=(""), ylab=("Cut Entitlements"), xlim=c(-0.25,1))
abline(v=0, col="grey", lty=3)
abline(v=mean(results$d.1+results$d.2), col="black", lty=1)
d<-density(results$d.1+results$d.2, kernel="gaussian")
polygon(c(0, rev(d$x[d$x<0])), c(0, rev(d$y[d$x<0])), col="maroon", border=NA)
plot(density(results$d.2, kernel="gaussian"), main=paste(""), xlab=(""), ylab=(""), xlim=c(-0.25,1))
abline(v=0, col="grey", lty=3)
abline(v=mean(results$d.2), col="black", lty=1)
d<-density(results$d.2, kernel="gaussian")
polygon(c(0, rev(d$x[d$x<0])), c(0, rev(d$y[d$x<0])), col="maroon", border=NA)
plot(density(results$d.3+results$d.2, kernel="gaussian"), main=paste(""), xlab=(""), ylab=(""), xlim=c(-0.25,1))
abline(v=0, col="grey", lty=3)
abline(v=mean(results$d.3+results$d.2), col="black", lty=1)
d<-density(results$d.3+results$d.2, kernel="gaussian")
polygon(c(0, rev(d$x[d$x<0])), c(0, rev(d$y[d$x<0])), col="maroon", border=NA)

plot(density(results$f.1+results$f.2, kernel="gaussian"), main=paste(""), xlab=(""), ylab=("Cut Defense"), xlim=c(-0.25,1))
abline(v=0, col="grey", lty=3)
abline(v=mean(results$f.1+results$f.2), col="black", lty=1)
plot(density(results$f.2, kernel="gaussian"), main=paste(""), xlab=(""), ylab=(""), xlim=c(-0.25,1))
abline(v=0, col="grey", lty=3)
abline(v=mean(results$f.2), col="black", lty=1)
plot(density(results$f.3+results$f.2, kernel="gaussian"), main=paste(""), xlab=(""), ylab=(""), xlim=c(-0.25,1))
abline(v=0, col="grey", lty=3)
abline(v=mean(results$f.3+results$f.2), col="black", lty=1)

plot(density(results$a.1+results$a.2, kernel="gaussian"), main=paste(""), xlab=(""), ylab=("Tax Rich"), xlim=c(-0.25,1))
abline(v=0, col="grey", lty=3)
abline(v=mean(results$a.1+results$a.2), col="black", lty=1)
plot(density(results$a.2, kernel="gaussian"), main=paste(""), xlab=(""), ylab=(""), xlim=c(-0.25,1))
abline(v=0, col="grey", lty=3)
abline(v=mean(results$a.2), col="black", lty=1)
plot(density(results$a.3+results$a.2, kernel="gaussian"), main=paste(""), xlab=(""), ylab=(""), xlim=c(-0.25,1))
abline(v=0, col="grey", lty=3)
abline(v=mean(results$a.3+results$a.2), col="black", lty=1)

plot(density(results$g.1+results$g.2, kernel="gaussian"), main=paste(""), xlab=(""), ylab=("Tax Both"), xlim=c(-0.25,1))
abline(v=0, col="grey", lty=3)
abline(v=mean(results$g.1+results$g.2), col="black", lty=1)
d<-density(results$g.1+results$g.2, kernel="gaussian")
polygon(c(0, rev(d$x[d$x<0])), c(0, rev(d$y[d$x<0])), col="maroon", border=NA)
plot(density(results$g.2, kernel="gaussian"), main=paste(""), xlab=(""), ylab=(""), xlim=c(-0.25,1))
abline(v=0, col="grey", lty=3)
abline(v=mean(results$g.2), col="black", lty=1)
d<-density(results$g.2, kernel="gaussian")
polygon(c(0, rev(d$x[d$x<0])), c(0, rev(d$y[d$x<0])), col="maroon", border=NA)
plot(density(results$g.3+results$g.2, kernel="gaussian"), main=paste(""), xlab=(""), ylab=(""), xlim=c(-0.25,1))
abline(v=0, col="grey", lty=3)
abline(v=mean(results$g.3+results$g.2), col="black", lty=1)
d<-density(results$g.3+results$g.2, kernel="gaussian")
polygon(c(0, rev(d$x[d$x<0])), c(0, rev(d$y[d$x<0])), col="maroon", border=NA)

plot(density(results$h.1+results$h.2, kernel="gaussian"), main=paste(""), xlab=("Correlation Parameter"), ylab=("Fed. Sales Tax"), xlim=c(-0.25,1))
abline(v=0, col="grey", lty=3)
abline(v=mean(results$h.1+results$h.2), col="black", lty=1)
plot(density(results$h.2, kernel="gaussian"), main=paste(""), xlab=("Correlation Parameter"), ylab=(""), xlim=c(-0.25,1))
abline(v=0, col="grey", lty=3)
abline(v=mean(results$h.2), col="black", lty=1)
plot(density(results$h.3+results$h.2, kernel="gaussian"), main=paste(""), xlab=("Correlation Parameter"), ylab=(""), xlim=c(-0.25,1))
abline(v=0, col="grey", lty=3)
abline(v=mean(results$h.3+results$h.2), col="black", lty=1)
title("Deliberative Persuasion", line=2, outer=TRUE, cex.main=2)

par(op)


resetGraph()





